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Abstract 

Hyperbolic- Relaxation  (HR)  systems  of  partial  differential  equations  (PDE’s) 
bear  the  promise  of  describing  rarefied  flows  in  the  transition  regime  (0.001  < 
Kn  <  10)  much  more  efficiently  than  pseudo-particle  methods  like  DSMC.  In 
order  to  arrive  at  a  maximally  efficient  numerical  methodology  a  Discontinuous- 
Galerkin  (DG)  spatial  discretization  was  combined  with  a  Hancock-type  tem¬ 
poral  integration  scheme.  A  linear  HR  system  of  two  equations,  describing 
1-D  advection-diffusion  in  the  equilibrium  limit,  was  used  as  a  test-bed  for 
error  analysis  and  numerical  experiments.  The  new  method  is  more  efficient 
than  any  of  the  traditional  DG  and  finite-volume  methods  tested,  and  re¬ 
mains  accurate  in  the  diffusion  limit.  The  method  is  ready  to  be  extended  to 
nonlinear  systems  and  multiple  dimensions.  In  addition  it  was  demonstrated 
that  spurious  inviscid  shocks  embedded  in  viscous  shock  structures,  found  on 
the  basis  of  HR  systems,  can  be  avoided  by  modifying  the  relaxation  term. 

Objective 

The  objective  of  the  current  research  project  is  to  develop  and  test  numer¬ 
ical  methods  for  Hyperbolic-Relaxation  (HR)  systems  of  PDE’s.  These  are 
first-order  systems  with  (usually  stiff)  source  terms  that  exhibit  hyperbolic 
behavior  (wave  propagation)  for  times  small  compared  to  the  relaxation  time 
r  (“frozen”  physics),  but  over  long  times  exhibit  diffusive  behavior  (“equilib¬ 
rium”  physics).  Such  systems  are  useful,  for  instance,  in  describing  the  flow 
of  rarefied  gases,  whether  in  the  upper  atmosphere  (re-entry,  braking)  or  in 
a  MEMS  device.  Numerical  methods  for  HR  systems  bear  the  promise  of 
describing  these  flows  much  more  efficiently  than  DSMC  (particle)  methods, 
for  Knudsen  numbers  as  high  as  10.  The  savings  could  be  orders  of  magni¬ 
tude  in  CPU  time  for  current  simulations,  greatly  reducing  the  turn-around 
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time  for  any  design  cycle. 

Crucial  to  the  success  of  the  HR  description  is  the  emergence  of  numeri¬ 
cal  methods  that  are  uniformly  accurate  with  regard  to  the  parameter  A t/r, 
where  At  is  the  time  step  used  in  the  scheme.  These  so-called  Asymptotic- 
Preserving  (AP)  methods  would  make  it  possible  to  compute  flows  with 
widely  varying  local  Knudsen  number. 

For  the  design  and  analysis  of  numerical  methods  we  use  an  HR  system 
of  only  two  equations,  the  so-called  “generalized  hyperbolic  heat  equation” 
(GHHE),  which  in  the  equilibrium  limit  reduces  to  a  single  advection-diffusion 
equation.  The  system  is  one-dimensional  and  is  taken  to  be  linear  when 
Fourier  analysis  is  in  order;  it  can  be  made  nonlinear  for  numerical  experi¬ 
mentation. 

Past  year’s  progress 

This  performance  report  covers  the  final  grant  period  9/1/2005  -  4/30/2006; 
it  also  serves  as  final  report. 

Results  on  DG  for  HR  systems 

We  wish  to  adopt  a  discretization  method  that  combines  compactness  with 
high-order  accuracy.  In  standard  finite-volume  methods,  higher-order  accu¬ 
racy  relies  on  piecewise-polynomial  reconstruction,  which  requires  extended 
stencils.  Discontinuous  Galerkin  (DG)  methods  overcome  the  draw-back  of 
reconstruction  by  using  extra  equations  for  updating  the  polynomial  repre¬ 
sentation  of  state  variables.  Currently,  the  most  successful  DG  methods  are 
semi-discretizations  combined  with  TVD  Runge-Kutta  (RK)  ODE  solvers, 
denoted  as  RKmDG(/c)  where  m  is  the  order  of  the  RK  method  and  k  is  the 
degree  of  the  polynomial  basis  functions. 

DG  methods  were  previously  shown  [1]  to  be  automatically  AP  provided 
the  advection  speed  is  small  compared  to  the  “frozen”  wave  speed.  This 
guarantee  is  lost  when  the  advection  speed  is  comparable  to  the  frozen  wave 
speed,  due  to  diffusive  numerical  errors  stemming  from  approximating  the 
advection  operator.  In  order  to  arrive  at  a  DG  method  that  remains  AP  for 
a  finite  advection  speed,  we  started  from  the  method  of  Arora  [2],  developed 
earlier  in  our  group,  which  uses  a  locally  implicit  characteristic  method  to 
evaluate  the  cell-interface  fluxes  with  strong  coupling  to  the  source  terms.  To 
realize  the  AP  property  for  all  practical  values  of  A f/r,  Arora  recommends 
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a  minimum  of  three  flux  evaluations  for  computing  the  flux  integral  over  one 
time  step.  Using  Huynh’s  [3]  idea  of  combining  a  Hancock-type  time  integra¬ 
tor  with  DG  spatial  discretization  (“upwind  moment  scheme”),  we  were  able 
to  cast  Arora’s  idea  in  a  form  that  uses  locally  implicit  updates  inside  a  cell, 
and  a  standard  Riemann  solver  to  couple  the  cells,  thus  greatly  simplifying 
the  method,  and  making  it  more  attractive  to  use. 

Our  HR  model  system  is  of  the  form 

dtu  +  <9xf(u)  =  “s(u),  (1) 

where  u  is  the  vector  of  conserved  variables,  f  is  the  flux  of  u,  s  is  the  source 
term,  and  e  is  the  nondimensional  relaxation  time.  The  numerical  results  we 
obtained  for  the  DG-Hancock  method  are  based  on  a  2  x  2  linear  HR  system 
with  u  =  [it,w]T,f  =  [v,u]T,  and  s  =  [0,ru  -  v]T  in  Eq.(l).  This  system 
has  “frozen”  wave  speeds  ±1  when  relaxation  is  weak  (e  »  1);  when  the 
relaxation  dominates  (e  <gC  1),  it  reduces  to  the  advection-diffusion  equation 
dtu  +  rdxu  ~  e(l  —  r2)dxxu,  with  an  “equilibrium”  wave  speed  of  r.  For 
stability,  |r|  <  1. 

The  solution  representation  is  piecewise  linear  (k  =  1);  thus,  the  gradient 
of  each  flow  variable  evolves  according  to  an  independent  update  equation. 
Details  of  the  method  and  the  numerical  experiments  are  presented  in  a 
conference  paper  by  Suzuki  and  Van  Leer  [4].  It  requires  solving  a  Riemann 
problem  twice  at  each  cell  interface  but  achieves  third-order  accuracy  in  time 
and  space.  We  solve  an  initial- value  problem  on  a  periodic  domain  with  a  har¬ 
monic  initial  condition  uq  =  vq  =  cos(27rx),  x  €  [0, 1];  the  other  parameters 
r,  e  are  chosen  so  that  the  reduced  equation  becomes  an  advection-dominated 
advection-diffusion  equation  (r  =  1/2,  e  =  10-5).  The  new  method  is  com¬ 
pared  with  a  second-order  Godunov-type  finite-volume  method,  HR2,  and  a 
DG(1)  method,  both  incorporating  the  IMEX-RK  method  [5].  At  t  =  300, 
when  the  wave  has  moved  150  times  its  own  length,  while  its  amplitude 
has  been  reduced  by  about  8%,  the  new  method  was  seen  to  be  the  least 
dissipative  and  dispersive  of  all,  whereas  the  RK2HR2  method  produces  a 
completely  inaccurate  solution.  A  grid-convergence  study  confirmed  that  the 
new  method  is  third-order  accurate,  as  expected  from  the  truncation-error 
analysis.  In  contrast,  RK2HR2  and  RK2DG(1)  show  second-order  conver¬ 
gence  in  the  L2-norm. 
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2.  Shock  representation  by  HR  systems 

Resolved  shock  structures  computed  on  the  basis  of  standard  hyperbolic- 
relaxation  systems  include,  for  most  Mach  numbers,  embedded  discontinu¬ 
ities  (inviscid  shocks)  that  are  not  validated  by  DSMC  calculations  [6].  The 
appearance  of  such  a  discontinuity  indicates  that  the  chosen  relaxation  term 
is  inadequate:  it  lacks  information  about  the  nonlinear  characteristic  fields 
of  the  hyperbolic  operator.  By  modifying  the  relaxation  term  we  have 
succeeded,  for  a  nonlinear  HR  model  system,  in  obtaining  the  correct  vis¬ 
cous  shock  profile.  In  order  to  show  that  a  time-marching  scheme  can  find 
the  proper  asymptotic  steady  solution  of  an  HR  system,  we  used  a  stan¬ 
dard  second-order  upwind-biased  finite- volume  scheme  with  a  two-stage  time- 
integrator  for  the  modified  system 


(2) 
(3) 

with  /  =  2.0,  r  =  0.1,  and  boundary  conditions 


ut  +  vx  =  0; 


,  ,2  2  (v-  |u2)/V 

vt  +  /Vu*  = - - - 


U2t 


u±OQ  =  ±U ,  u±00  =  U,  U=  1.  (4) 

The  original  relaxation  term  was  (v  —  \u2)/t.  Figure  1  shows  the  numeri¬ 
cal  results  plotted  on  top  of  the  exact  profile,  for  Ax  =  0.075.  There  is  no 
trace  of  an  inviscid  jump,  and  the  agreement  appears  to  be  excellent.  For 
comparison,  Figure  2  shows  the  incorrect  profile  obtained  with  the  original 
system;  it  has  an  infinite  derivative  at  the  origin).  A  grid-refinement  study 
confirms  the  second-order  convergence  of  the  numerical  to  the  exact  solution. 
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Figure  1:  Steady  Burgers  shock 
profile  (line,  exact  solution  cell 
averaged)  and  numerical  approx¬ 
imation  (symbols)  obtained  with 
the  modified  hyperbolic-relaxation 
system  (3);  r  =  0.1,  Ax  =  0.075. 
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Figure  2:  Steady  Burgers  shock 
profile  and  numerical  approxima¬ 
tion  obtained  with  the  original 
hyperbolic- relaxation  system;  r  = 
0.1,  Ax  =  0.075.  The  numerical 
profile  is  too  steep;  its  derviative 
in  the  origin  is  infinite. 
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